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Introduction 


This  report  details  work  on  applying  real  gas  properties  to  the  determination  of  critical  values  of 
shock  wave  angle  for  regular  reflections  (aN )  and  Mach  reflections  (ad  ).  An  excellent  review  of 
regular  and  Mach  reflections  has  been  given  in  Hornung  [1],  which  describes  the  physical  process. 
More  recent  work  has  been  done  with  respect  to  numerical  simulations  in  the  dual-solution-domain  [2], 
where  either  regular  or  Mach  reflections  can  exist.  Most  of  the  work  to  date  has,  however,  focused  on 
relatively  weak  shock  waves  and  perfect  gas,  constant-  y  shock  wave  processes.  The  purpose  of  the 
current  research  is  to  extend  these  results  to  Mach  numbers  up  to  18,  where  variable  properties  and 
reactions  can  be  important. 

There  are  two  main  effects  when  considering  high  Mach  number  flows  with  strong  shock 
waves.  The  first  effect  of  importance  is  a  variable  specific  heat.  As  higher  temperatures  are  achieved 
behind  the  shock,  the  specific  heat  for  molecules  will  typically  increase  as  vibrational  and  eventually 
electronic  states  become  excited  in  the  molecules,  shifting  the  specific  heat  ratio  towards  one  and 
effecting  the  shock  wave.  The  second  main  effect  is  reactions.  Reactions  can  occur  behind  the  shock 
due  to  combustion  or  dissociation,  which  either  add  or  remove  heat  from  the  gas  behind  the  shock 
wave.  For  this  report  we  focus  on  dissociation,  and  use  an  equilibrium  calculation  to  determine  the 
shocked  gas  composition. 


Problem  Statement 

The  problem  of  interest  is  illustrated  in  Figure  1.  Flow  enters  the  domain  at  a  given 
temperature,  pressure,  mole-fractions,  and  velocity.  For  the  cases  described  here,  the  inflow  pressure 
and  temperature  are  one  atm  and  298.15  K.  The  flow  encounters  a  wedge,  which  creates  an  oblique 
shock  wave  (we  only  consider  angles  producing  oblique,  attached  shock  waves).  This  shock  wave  is 
then  reflected  off  of  the  top  wall  or  symmetry  plane.  Two  types  of  reflections  can  occur,  dependent  on 
the  inflow  conditions  and  the  wedge  angle.  The  first  type  of  reflection  is  a  regular  reflection,  where  an 
oblique  shock  wave  is  reflected  off  of  the  surface.  The  second  type  is  a  Mach  reflection,  where  the 
incident  shock  transitions  to  a  normal  shock  near  the  top  wall  or  reflected  boundary.  For  inflow  Mach 
numbers  greater  than  2.2,  there  are  two  critical  values  for  the  wedge  angle  (and  thus  shock  wave  angle) 
that  determine  the  type  of  reflection  that  occurs.  For  a<aN  only  regular  reflections  occur,  and  for 
a  >  ad  only  Mach  reflections  occur.  However,  for  aN  <  a  <  ad  both  regular  and  Mach  reflections 
occur.  This  region  is  known  as  the  dual-solution-domain  (DSD).  Both  critical  values  for  a  are  found 
by  considering  the  regular  reflection  case. 


Figure  1.  Schematic  of  regular  reflection  (left)  and  Mach  reflection  (right)  off  of  a  solid  wall  boundary. 
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For  regular  reflections,  the  analysis  is  divided  into  two  oblique  shock  wave  solutions. 
Conditions  ahead  of  the  incident  shock  wave  are  described  as  condition  (1),  and  behind  the  incident 
shock  wave  we  have  condition  (2).  These  conditions  serve  as  the  pre-shock  conditions  for  the  reflected 
shock  wave  calculation.  The  flow  after  the  reflected  shock  wave  is  represented  as  condition  (3).  The 
next  section  details  the  solution  method  used  for  oblique  shock  waves. 


Oblique  Shock  Wave  Analysis 


The  oblique  shock-wave  solution  is  found  by  applying  the  appropriate  mass,  momentum,  and 
energy  conservation  equations  along  with  an  equation  of  state  and  thermodynamic  relations.  For  an 
oblique  shock- wave,  we  first  identify  the  component  of  velocity  normal  to  the  shock- wave  ( u)  and  the 
component  of  velocity  tangential  to  the  shock-wave  ( w).  Labeling  pre-shocked  properties  with  a 
subscript  1  and  shocked  properties  with  a  subscript  2,  we  can  develop  conservation  equations  for  the 
shocked  gases  after  an  oblique  shock-wave.  For  a  shock  wave  angle  of  a,  the  velocity  components  are 

=  Vasina  (1) 

Wj  =  VjCosa  (2) 

where  Vj  is  the  initial  velocity  of  the  flow.  Conservation  of  mass,  momentum,  and  energy  reduces  to: 


=  p2U2 

(3) 

<N 

* 

II 

ST 

(4) 

Pi  +  PlUl  ~  P2  +  PlU2 

(5) 

2  2 

,  Mi  Up 

h]  H - —  h1  H - 

(6) 

2  2 

In  addition,  we  have  auxiliary  equations  in  the  form  of  the  equation  of  state  and  thermodynamic 
relations  to  close  the  system  of  equations.  This  can  be  summarized  as: 

Pi  =  P^PiJ^i)  (7) 

After  computing  the  shocked  gas  properties,  we  can  calculate  the  deflection  angle  as 

6  =  a-  tan-1—  (8) 


where  6  is  the  deflection  angle.  For  a  perfect  gas  with  a  constant  y,  the  above  relations  can  be  reduced 
to  two  independent  variables,  Mx  and  a.  Pressure  and  temperature  ratio,  and  Mach  number  are  all 
expressed  as  functions  of  Mx  and  y.  For  real  gases  with  possible  reactions,  this  is  no  longer  possible, 
and  the  resultant  pressure  and  temperature  ratios  are  dependent  on  the  pre-shocked  pressure, 
temperature,  velocity,  and  mixture  composition.  To  solve  this  system,  we  need  to  know  enthalpy  and 
specific  heat  as  a  function  of  temperature  and  pressure.  For  the  equilibrium  reaction  computations,  we 
also  need  to  determine  mixture  composition  as  a  function  of  pressure  and  enthalpy.  Appendix  A  gives 
more  details  on  the  iterative  solution  to  the  conservation  equations,  Appendix  B  describes  how  the 
thermal  properties  of  the  gas  are  computed,  and  Appendix  C  describes  how  the  equilibrium 
composition  of  the  gas  is  found. 
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Deflection  Diagrams 


The  preceding  equations  represent  the  core  computations  for  our  problem  of  interest.  For  our 
problem,  we  are  interested  in  understanding  how  the  pressure  ratio  ( p/p1 )  and  deflection  angle  ( 9 ) 
vary  as  the  incident  angle  varies  (a).  As  mentioned  previously,  for  a  perfect  gas  with  no  reactions, 
these  variables  are  dependent  solely  on  the  upstream  Mach  number.  However,  for  variable  specific 
heats  and  reactions,  this  problem  is  dependent  on  all  of  the  inflow  conditions,  u1,PvT1,Xl  i  =  C(1).  In 
the  computations  completed  in  this  report,  the  upstream  pressure  and  temperature  are  one  atm  and 
298.15  K,  and  the  mole-fractions  are  specified  as  either  air  or  detonation  products.  We  still  report  the 
findings  in  terms  of  the  upstream  Mach  number,  to  provide  comparisons  with  previous  work  with 
constant  specific  heat  ratios.  Note,  however,  that  if  the  inflow  temperature,  pressure,  or  mole-fractions 
are  altered,  the  curves  presented  will  also  shift. 

The  minimum  incident  angle  that  will  produce  a  shock- wave  is  based  on  the  requirement  that 
the  velocity  normal  to  the  shock- wave  must  at  least  be  sonic.  This  reduces  to  the  equation: 


a.  =  sin 

min 


-1 


M, 


(9) 


a  is  then  varied  between  amin  and  jt /2,  which  represents  a  normal  shock.  Solving  the  oblique  shock 
relations  as  a  varies  between  amin  and  Jt/2  gives  us  deflection-pressure  diagrams.  In  Figure  2,  we 
consider  inflow  air  at  standard  temperature  and  pressure  and  a  velocity  corresponding  to  Mach  3.  The 
independent  variable  to  produce  this  plot  is  the  incident  angle,  a,  and  the  reflected  angle,  /3,  even 
though  the  horizontal  axis  for  these  plots  is  the  deflection  angle,  6. 

An  important  characteristic  of  pressure-deflection  diagrams  is  that  we  can  represent  the  entire 
incident-reflected  shock- wave  phenomenon  on  one  plot.  We  simply  superimpose  both  pressure- 
deflection  diagrams  on  the  same  plot,  translating  the  reflected  shock-wave  diagram  to  the  correct  initial 
pressure  and  deflection  angle.  The  pressure-deflection  diagram  in  Figure  2  represents  the  shock- wave 
pattern  illustrated  in  Figure  1.  Conditions  at  location  (3)  are  found  by  determining  where  the  pressure- 
deflection  diagram  for  the  reflected  shock- wave  intersects  the  vertical  axis.  Note  that  each  deflection 
angle  represents  two  unique  values  for  pressure  and  also  for  shock  wave  angle,  a,  and  that  for  any 
given  Mach  number,  there  is  a  maximum  deflection  angle  for  which  an  attached  oblique  shock- wave 
can  exist.  More  details  for  deflection-pressure 
diagrams  is  given  in  the  review  article  [1].  2.5 

In  a  like  manner,  we  can  show  a  series  of 
deflection-pressure  diagrams  for  the  same  inflow  2 
conditions  (1),  but  differing  deflection  angles  6.  ^ 

Q_  -|  cj 

This  is  shown  in  Figure  3  for  two  different  Mach  5; 
numbers,  Mach  3  and  8.  There  are  a  some  — 

important  things  to  note  in  Figure  3.  First,  as  the 
deflection  angle  increases  (and  thus  the  shock  0  5 

wave  angle  a  increases),  the  incident  shock  wave 
becomes  stronger,  and  the  reflected  shock  wave  0 

becomes  weaker  due  to  the  lower  Mach  number  in 
region  (2).  This  results  in  a  smaller  deflection- 
pressure  curve  for  the  reflected  shock  wave,  as  can  Figure  2.  Pressure-deflection  diagram  for  incident  and 

be  seen  in  the  Figure  There  are  two  critical  values  reflecte<^  shock.  Upstream  condition  is  air  at  standard 

.  ......  ...  temperature  and  pressure,  and  the  Mach  number  of  the  flow  is  3 

for  a  that  are  labeled  in  the  Figure.  The  first  and  y= 1.4. 

critical  value,  labeled  aN ,  is  where  the  pressure  at 
location  (3)  equals  the  normal  shock  pressure  for 
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Figure  3.  Deflection-pressure  diagrams  for  incident  and  reflected  shock-waves.  Inflow  Mach  number  3  (left)  and  8  (right), 
y=1.4.  Pressure -deflection  diagrams  for  critical  values  of  the  incident  shock-wave  angle  are  given  by  dashed  lines. 


the  inflow  conditions.  This  normal  shock  pressure  occurs  at  a=  jz  / 2  and  9  =  0  on  the  incident  shock- 
wave  deflection  diagram.  For  values  of  a  less  than  aN ,  a  regular  reflection  will  always  occur.  The 
second  critical  value,  labeled  ad  ,  is  where  the  maximum  deflection  angle  of  the  reflected  shock- wave 
is  exactly  equal  to  the  deflection  angle  of  the  first  shock-wave.  For  a  greater  than  ad  ,  a  Mach 
reflection  will  always  occur.  The  region  between  the  two  critical  values,  where  aN  <  a  <  ad ,  is  known 
as  the  dual-solution-domain  (DSD),  either  a  regular  reflection  or  a  Mach  reflection  can  exist  in  the 
DSD.  In  general,  as  the  Mach  number  is  increased,  the  difference  between  aN  and  ad  also  increases. 

The  purpose  of  this  report  is  to  investigate  how  these  critical  values  for  a  vary  at  high  Mach 
numbers  where  variable  thermal  properties  and  dissociation  become  important.  For  perfect  gases,  the 
critical  a's  are  only  a  function  of  upstream  Mach  number  and  specific  heat  ratio.  This  allows  the 
expressions  for  determining  these  values  to  be  simplified  dramatically.  However,  because  we  do  not 
make  the  assumption  of  perfects  gases,  we  need  to  use  iterative  techniques  in  order  to  determine  the 
correct  values  for  aN  and  ad  ,  using  the  criteria  mentioned  previously  to  determine  the  critical  values. 
For  aN ,  this  is: 

^(aiV’C(i))  -  -Pn(Ca))  =  0  (10) 

where  Pf  x,C(l  )  is  the  pressure  at  station  (3)  in  Figure  1,  as  a  function  of  the  incident  shock-wave 
angle  and  upstream  conditions  represented  by  C(1) .  F„(C(1))  is  the  normal  shock- wave  pressure  jump 

using  the  conditions  at  station  (1)  upstream  of  the  incident  shock- wave.  The  Mach  reflection  limit, 
ad  ,  is  solved  by 

^(ad’C(i))  “  0m<„(C(2))  =  0  .  (11) 

The  first  set  of  solutions  for  aN  and  ad  are  shown  in  Figure  4.  This  solution  was  done  for  a 
perfect  gas  with  y  =  1.4,  using  the  method  developed  for  gases  with  variable  specific  heats  and  possible 
reactions.  The  results  are  compared  with  Khotyanovsky  et  al  [2],  and  shows  that  the  methods  for 
determining  aN  and  ad  can  reproduce  previous  results  for  perfect  gases.  Since  our  interest  is  in  the 
high  Mach  number  range,  the  Mach  number  range  below  3  was  not  investigated  for  this  report.  Their 
results  were  limited  to  Mach  numbers  below  6.5.  Results  for  this  paper  examine  critical  incident 
shock- wave  angles  for  Mach  numbers  in  the  range  of  3.0  all  the  way  to  18.0,  to  correspond  to  the  large 
range  of  conditions  seen  with  blast  waves. 
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Figure  5  shows  results  for  aN  and  ad 
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for  perfect  gases,  where  y  varies  from  1.2  to  "5. 

1.6.  For  comparison,  we  also  compute  the  curve  < 
for  curve-fit  air  properties  on  the  same  plot.  £ 

Although  for  such  a  wide  range  of  Mach  g 

numbers,  we  expect  temperatures  to  be  high 
enough  to  require  variable  properties  for 
accuracy,  the  trends  that  we  see  for  the  perfect 
gas  cases  are  instructive.  In  this  plot,  we  see 
that  as  y  increases,  aN  tends  to  increase  and  ac 
tends  to  decrease,  creating  a  smaller  dual¬ 
solution-domain  region.  The  value  for  y  has  a 
larger  effect  on  the  value  of  ad  than  aN  .  To 
relate  this  to  real  gases,  note  that  as  the 
temperature  rises,  the  value  for  Cp  generally 
increases,  which  will  shift  the  value  towards  a 
lower/.  This  is  seen  in  Figure  5  for  air  with 
variable  properties,  which  has  a  value  of  1.4  at  low  temperatures,  but  will  increase  as  the  temperature 
rises.  For  detonation  and  combustion  products,  the  value  for  /  tends  to  be  lower  than  1.4  even  at  low 
temperatures,  due  to  the  larger  heat  capacity  of  combustion  products.  Therefore  we  expect  the 
detonation  results  to  follow  the  lower  /  curves  rather  than  the  higher  /  curves.  The  final  interesting 
result  from  this  figure  is  how  the  variable  properties  curves  vary  from  the  perfect  gas  /  =  1.4  curves. 
There  is  very  little  difference  for  aN ,  even  in  the  high  Mach  number  limit.  There  is,  however,  more 
variation  for  the  ad  values.  This  is  not  surprising,  since  the  larger  incident  angles  result  in  much  larger 
temperature  jumps  across  the  shock,  and  would  thus  have  a  comparatively  larger  change  in  /. 

The  second  main  result,  shown  in  Figure  6,  examines  the  critical  values  for  incident  shock- wave 
angle  for  different  combustion  type  models.  For  comparison,  we  have  placed  the  variable  properties, 
air  curves  on  this  figure  as  well.  For  the  detonation  products,  we  use  an  approximation  to  the  products 
found  in  Volk  and  Schedlbauer  [3]. 


Figure  4.  Critical  values  for  incident  shock-wave  angle  for  a 
perfect  gas,  y=1.4.  Symbols  represent  data  from  Khotyanovsky  et 
al[2 ]. 


Species 

Mole-fraction 

Specific  heat  ratio  at  300  K 

Carbon  Monoxide  (CO) 

0.174 

1.400 

Carbon  Dioxide  (C02) 

0.100 

1.289 

Nitrogen  (N2) 

0.136 

1.400 

Steam  (H20) 

0.198 

1.327 

Hydrogen  (H2) 

0.034 

1.409 

Carbon  graphite  (C) 

0.359 

1.000 

Each  of  the  species  is  curve-fit  from  200  K  to  6000  K.  Above  6000  K,  we  could  not  find  curve-fits  for 
some  of  the  species  (H20  and  C),  and  instead  extended  the  specific  heat  as  a  constant  above  6000  K. 
This  will  effect  only  the  high  Mach  number  range  of  ad  calculation,  about  Mach  13.2.  Below  this,  and 
for  all  of  the  values  for  aN ,  the  temperature  is  within  the  6000  K  range. 
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Figure  5.  Critical  values  for  incident  shock-wave  angle  for 
perfect  gases  and  air  with  variable  properties.  aN  is 
represented  as  a  solid  line,  and  ad  as  a  dashed  line. 


Figure  6.  Critical  values  for  incident  shock-wave  angle  for  air 
and  detonation  products  using  variable  properties.  For 
comparison,  two  perfect  gas  cases  are  also  presented. 


With  this  formulation,  we  also  have  the  ability  to  examine  three  different  perfect  gases  at 
stations  (1),  (2),  and  (3)  in  Figure  1,  as  well  as  temperature-dependent  variable  properties.  Thus,  at 
each  station,  we  can  have  a  different  value  for  y.  This  could  either  represent  an  increase  in  specific 
heats  due  to  higher  temperatures  after  the  incident  shock  wave,  or  it  could  also  represent  different 
compositions  due  to  chemical  reactions.  In  this  case,  we  look  at  a  case  where  the  y  varies  from  1.4  in 
front  of  the  incident  shock  wave,  to  1.2  behind  both  the  incident  and  reflected  shock  waves.  In  Figure 
6,  we  see  that  the  value  for  aN  appears  to  be  largely  dependent  on  the  specific  heat  ratio  after  the 
incident  shock  wave.  The  perfect  gas  solutions  for  y  =  1.2  after  the  incident  shock  wave  both  do  a 
remarkable  job  of  reproducing  the  variable  properties  result  with  combustion  products.  Again,  this  is 
most  likely  due  to  the  minimal  temperature  rise  after  shock  waves  with  such  low  angles  of  incidence, 
even  at  high  Mach  numbers.  More  variation  is  seen  with  the  ad  computation,  especially  with  the 
combustion  products.  Again,  this  is  a  similar  trend  as  what  was  seen  in  Figure  5  with  variable 
properties  compared  with  y  =  1.4,  likely  due  to  increasing  specific  heats  as  the  temperature  increases 
for  the  variable  property  calculations.  In  all  cases  seen  in  Figure  6,  the  y  =  1.4  perfect  gas  solution  is 
considerably  different  than  any  of  the  other  approximations,  and  shows  that  it  is  important  to  consider 
the  effect  of  different  y’s  and  variable  properties. 

The  final  result  that  we  present  is  shown  in  Figures  7  and  8.  In  these  plots,  we  compare  a 
perfect  gas  with  y  =  1.4,  frozen  (non-reacting)  air  with  variable  properties,  and  air  with  equilibrium 
reactions  and  variable  properties.  We  show  both  the  critical  shock  wave  angles  as  Mach  number 
increases,  and  also  the  temperature  found  behind  the  reflected  shock  wave,  which  represents  the  highest 
temperature  in  the  system.  For  the  air  with  equilibrium  reactions,  we  consider  the  species  N2,  02,  N, 

O,  and  NO  as  the  main  constituents.  Results  show  again  that  aN  has  only  minor  variations  when  we 
consider  variable  properties  and  equilibrium  reactions,  but  ad  has  a  much  more  pronounced  variation. 
Also  note  that  the  value  for  ad  increases  as  we  consider  equilibrium  reactions.  Intuitively,  we  might 
think  that  the  value  should  actually  decrease,  since  the  lower  temperatures  and  more  monatomic  nature 
of  the  gas  suggests  a  higher  value  for  y.  This,  however,  is  countered  by  the  dissociation  reactions 
extracting  heat  from  the  shocked  gases,  which  tends  to  shift  the  value  for  ad  upwards.  Examining 
Figure  8,  we  can  see  the  significance  of  the  heat  extraction,  as  the  temperature  is  reduced  from  12,674 
K  to  7458  K  due  to  the  dissociation  of  the  molecules.  Since  dissociation  had  only  a  minor  effect  on  the 
value  for  aN  with  air,  we  have  decided  not  to  do  equilibrium  reactions  for  the  detonation  products. 
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Figure  7.  Critical  values  for  incident  shock  wave  angle  for 
air  with  variable  properties ,  and  with  and  without 
dissociation. 


Figure  8.  Temperature  after  the  reflected  shock  at  the  critical 
values  ad  (dashed)  and  aN  (solid). 


This  would  prove  much  more  difficult  due  to  the  range  of  possible  species  and  the  presence  of  solid 
graphite  within  the  flow. 


Conclusions 

This  short  report  has  described  calculations  performed  to  determine  the  critical  values  for  the 
shock  wave  angle  aN  and  ad  for  Mach  numbers  up  to  18,  for  both  air  and  detonation  products.  We 
have  found  that  both  critical  values  are  very  sensitive  to  the  specific  heat  ratio  after  the  incident  shock 
wave.  For  perfect  gases,  as  y  decreases,  the  value  for  ad  increases  and  the  value  for  aN  decreases, 
thus  increasing  significantly  the  dual-solution-domain  where  both  regular  and  Mach  reflections  can 
exist.  For  aN ,  the  results  were  fairly  insensitive  to  variable  specific  heats  even  at  very  high  Mach 
number,  and  adequate  approximations  can  be  made  without  considering  a  variable  specific  heat, 
although  careful  attention  should  be  made  to  the  selection  of  y.  ad  showed  a  much  greater  sensitivity 
to  variable  properties,  and  both  variable  properties  and  equilibrium  reactions  should  be  considered 
when  computing  this  critical  value  even  at  relatively  low  Mach  numbers.  Although  the  value  for  aN  is 
fairly  insensitive  to  variable  properties,  the  temperature  behind  the  reflected  shock  does  vary 
significantly  for  variable  specific  heats.  If  this  parameter  is  important  for  the  specific  application  under 
investigation  than  variable  specific  heats  should  be  considered.  Equilibrium  chemistry  was  not 
included  with  the  detonation  products  because  our  main  interest  was  determining  good  values  for  aN  at 
very  high  Mach  numbers.  The  complete  iterative  equations,  as  well  as  the  curve-fits  used  for  the 
specific  heat,  enthalpy,  and  entropy  computations,  are  given  in  the  Appendix. 
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Appendix  A.  Solution  Procedure  for  Real  Gases 

To  compute  properties  across  an  oblique  shock  for  variable  specific  heat  properties,  we  must 
use  an  iterative  algorithm.  The  algorithm  is  based  on  the  conserved  variables  across  a  shock.  For  an 
oblique  shock  where  the  normal  velocity  across  the  shock  is  and  the  tangential  component  is  w,  the 
mass,  momentum,  and  energy  equations  can  be  reduced  to: 

We  can  rearrange  these  equations  in  terms  of  an  iterative  variable  R  =  pj  p2.  This  variable  will 
always  be  less  than  1.  The  momentum,  energy,  and  continuity  equations  can  be  rewritten: 


Pl  =  Pl+P,ul{l-R) 

(12) 

h  =  K  +  ^-R2) 

(13) 

u2  =  u{R 

(14) 

The  iterative  procedure  employed  is  relatively  simple.  First,  an  initial  value  for  R  is  guessed. 
There  are  several  ways  to  get  a  good  guess  for  R,  but  what  we've  found  is  robust  convergence  is  fairly 
independent  of  this  initial  guess.  Therefore,  we  typically  assume  an  initial  value  of  R  =  0.1.  The 
iterative  process  is  then: 

1.  Solve  equation  (12)  for  p2 . 

2.  Solve  equation  (13)  for  h2. 

3.  Calculate  p2  from  equation  (7). 

4.  Calculate  new  Rn+i  =  pj p2 . 

5.  Check  for  convergence  of  R'*' .  If  not  converged,  go  back  to  step  1. 

6.  Calculate  u2  from  equation  (14). 

We  converge  -  /?")  //?"*"'  down  to  1  x  10"8 . 

For  non-reacting  shock  waves  (where  no  dissociation  occurs),  equation  (7)  is  fairly  simple  to 
compute.  For  reacting  shock  waves,  the  situation  can  be  more  complex.  If  we  know  the  composition 
of  the  products,  such  as  assuming  a  fuel  gas  is  completely  consumed  to  products,  then  we  use  a  similar 
iterative  solution  as  above.  If,  however,  we  are  interested  in  equilibrium  concentrations  (for  instance,  if 
we  are  assuming  dissociation  of  air  species),  then  we  must  solve  an  iterative  set  of  equations  for  the 
composition  and  temperature,  holding  enthalpy  and  pressure  constant. 

In  both  cases,  once  temperature  and  compositions  are  found,  we  use  the  ideal  gas  law  to 
compute  the  density: 

NS  NG 

P  =  jf2x.M</2x>  (15) 

(=1  (=1 

where  NS  is  the  total  number  of  species,  and  NG  is  the  number  of  gas  species.  If  all  species  are 
gaseous,  than  =1. 
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Appendix  B.  Thermal  Properties  Calculation 

Thermal  properties  for  the  gas  mixture  is  obtained  through  curve-fits  for  the  enthalpy,  entropy, 
and  specific  heat.  These  properties  take  the  general  form  for  species  j: 

c a*) 

For  the  CEA  database  [4],  each  species  specific  heat  is  modeled  by  three  7th-order  polynomials  in  the 
form  of  Eqn.  16,  with  temperature  ranges  from  200-1000  K,  from  1000-6000  K,  and  6000-20,000  K. 
The  value  for  mi  varies  from  -2  to  5.  Enthalpy  and  entropy  values  are  obtained  by  using  the  relations: 

h°(T)  =  fCpj(r)dT+cjN+l  (17) 

s°(T)  =  fC±J^-dT+ci„tl  (18) 

For  each  curve-fit,  enthalpy  and  entropy  at  any  temperature  is  specified  with  the  above  definitions  and 
a  constant  representing  the  enthalpy  and  entropy  of  formation.  For  a  mixture,  the  specific  enthalpy  and 
entropy  (molar-basis)  are  simply  summations  of  species  enthalpies  multiplied  by  their  mole-fractions, 
i.e., 

h(T)  =  2x,K‘(T)  (19) 

For  shock-wave  problems,  we  typically  know  the  specific  enthalpy  and  composition,  and  need  to  solve 
for  temperature  T .  In  these  cases,  we  use  an  iterative  Newton-Rhapson  method  to  converge  on  the 
appropriate  temperature,  which  works  well  due  to  the  smoothness  of  the  specific  heat  curves.  The  bar 
above  the  variables  indicate  that  these  are  mole  based  properties;  to  convert  to  mass  based  properties 
(that  are  used  in  the  conservation  equations)  we  divide  by  the  molecular  weight.  For  species  in  the 

CEA  database,  the  C°Pj  is  fitted  using  a  7th-order  polynomial,  as  shown  below: 


C  j  l  Cj  2 

jr  +  ^r+cl3  +  cp 


,£al 

T 


i(T)  =  -^l~C-f  +  cj,3\nT. 


Il 

T 


r4 

(20) 

T 5 

•7T+C;’8 

(21) 

T4 

(22) 

- *" C  i  9 

4 
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Appendix  C.  Chemical  Equilibrium 

The  equilibrium  calculation  is  discussed  in  detail  in  NASA  Reference  Publication  1311  [4],  We 
only  summarize  the  equations  here  for  completeness.  We  do  not  use  the  CEA  program  directly,  but 
have  instead  re-implemented  the  algorithm  to  better  incorporate  it  into  the  current  solution  procedure. 
The  equilibrium  calculation  minimizes  the  Gibbs  free  energy  for  the  system.  This  minimization  is  done 
by  solving  the  system  8G  =  0,  and  using  a  Newton-Rhapson  iterative  scheme  to  converge  the  system. 
By  applying  the  constraint  that  the  element  number  must  be  conserved,  we  get  a  set  of  reduced  Gibbs 
iterative  equations  assuming  no  condensed  phases: 


NE  NG 

22 

i= 1  7=1 


(  NG  \  ( NG  „  _  ko\ 


a ,  a  n  -it-  + 

y  j  i 
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\  7-1 
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RT 
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NE  is  the  number  of  elements,  nj  is  defined  as  the  number  of  moles  of  species  j  per  kg  mixture.  atj 
are  the  number  of  atoms  of  element  i  in  species  j.  bj  =  ^  aijnj .  bi  o  indicates  the  assigned  number  of 
moles  of  atom  i  per  kg  mixture,  n  is  the  moles  of  gas  per  kg  of  mixture,  ^  nj .  uj  is  the  chemical 
potential  of  species  j,  defined  below  for  gases: 

Uj  =  hj  -s°T  +  RTln^n j  / nj  +  RTln^p / p0')  (26) 

where  p°  is  the  standard- state  pressure.  ni  is  related  the  Lagrange  multiplier  that  is  needed  for  the 
solution  method,  where  Jti  =  -XJRT . 

The  iteration  variables  for  Eqns  (23-25)  are  nx,  Ain n ,  and  AlnT.  This  gives  us  NE  +  2 
equations  to  solve  simultaneously.  We  use  the  LAPACK  routines  to  solve  this  system  of  equations. 
The  gaseous  species  concentrations  are  then: 

„  NE  foo 

Ain n = - —  +  Y  a.  jt,  +  Alnn  -l — —  AlnT  (27) 

RT  RT 

We  then  use  these  values  to  iterate  on  the  composition  and  temperature: 

lnn^  =  lnrt(;,)  +  A(m)(Alnn  .)(m)  (28) 

lnn(m+1)  =  lnn(m)  +  A(m)(Aln  njm)  (29) 

lnT(m+1)  =  lnT(m)  +  A(m)(AlnT)(m)  (30) 

where  A('")  is  a  relaxation  factor  to  ensure  that  the  corrections  remain  constrained  for  the  mth  iteration. 
Reference  4  gives  a  detailed  description  of  how  A*"'-1  is  chosen.  The  convergence  criteria  is: 
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where  ntof  =  ^ tij  and  refers  to  the  chemical  element  with  the  largest  value  of  bio. 

An  initial  guess  for  the  composition  is  provided  by  the  user,  along  with  an  initial  guess  for  the 
temperature.  For  the  shock  wave  calculations  presented  in  this  report,  we  assume  no  dissociation  to 
obtain  an  initial  guess  for  the  composition  and  temperature.  Convergence  usually  occurs  within  18  to 
19  iterations  for  the  air  equilibrium  species.  This  convergence  is  clearly  dependent  on  the  problem, 
however,  and  the  number  of  species  and  elements  considered  in  the  system. 
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